########################################################################################################################
## This program calculates the inference of estimates and counterfactuals ##
## Model refers to "Distribution Regression with Selection"
## Authors: Chernozhukov, V., I. Fernandez-Val and S. Luo
## Last Update: 6/17/2018

## 80s vs 90s vs 100s within gender 

## male vs female within year(s)




rm(list=ls());

library(parallel);
# library(sampleSelection);
# library(maxLik);
library(mvtnorm);
library(abind);
library(MASS);
library(pryr);

#####################
# Extract Arguments #
#####################
#Define default values for variables
n.boot         <- 500;
counter.method <- 2; 
# 1 for counterfactual decades within gender; 
# 2 for counterfactual sex within time period; 
# 3 for counterfactual 18 years within gender;
rho.spec       <- 1; 
# Specification of rho: 1 for default; 
#                       2 for education-specific rho's; 
#                       3 for couple-specific rho's;
#                       4 for time-specific rho's;
#                       5 for linear expression of X;
#                       6 for age; 
#                       7 for number of kids;
#                       8 for linear expression of years of schooling;
#                       9 for 4 education dummies;
#                       10 for 4 education dummies x couple;
#                       11 for age x couple;
#                       12 for 4 age dummies;
#                       13 for numkids x couple;
#                       14 for 4 age dummies x couple;
#                       15 for linear in t;
#                       16 for quardratic in t;
#                       17 for yrseduc x couple;
#                       18 for linear in t x couple;
#                       19 for 6 time dummies;
#                       20 for propensity score of employment;
#                       21 for pscore x couple;
#                       22 for beninc0 + m_inc0;
sample.id      <- 11378;
l.trim         <- 10;
u.trim         <- 90;
alpha          <- 0.05;

#Get aruments from the command line
argv <- commandArgs(TRUE);

# Check if the command line is not empty and convert values to numerical values
if(length(argv) > 0){
  n.boot         <- as.numeric(argv[1]);
  counter.method <- as.numeric(argv[2]);
  rho.spec       <- as.numeric(argv[3]);
  sample.id      <- as.numeric(argv[4]);
  l.trim         <- as.numeric(argv[5]);
  u.trim         <- as.numeric(argv[6]);
  alpha          <- as.numeric(argv[7]);
};


str.sample <- ifelse(counter.method==1, ifelse(sample.id==1, "_m","_f"), 
                     ifelse(counter.method==2, paste("_",sample.id,sep=""), ifelse(sample.id==1, "_mhalf","_fhalf")));
str.rho.all <- c("","_educ","_couple","_time","_allx","_age","_numkids","_yrseduc","_4dum","_4dumxcouple",
                 "_agexcouple","_4agedum","_numkidsxcouple","_4agedumxcouple","_lintime","_quartime","_yrseducxcouple",
                 "_lintimexcouple","_6tdum","_pscore","_pscorexcouple","_iv");
str.rho <- str.rho.all[rho.spec];

str.trim <- paste("_", ifelse(l.trim<10, "0",""), l.trim, ifelse(u.trim<10, "0",""), u.trim, sep="");

str.cb <- ifelse(alpha==0.05, "","_c90");

str.robust <- ""; #"_blundell3";

str.boot <- ifelse(n.boot==200,"",paste("_",n.boot,sep=""));


###########################
# Load estimation results #
###########################
load(paste("Est",str.sample,str.rho,str.robust,".RData",sep=""));
source("DRS_Source.R");


####################
# Parallel Setting #
####################

# Calculate the namber of cores
n.cores <- ceiling(as.numeric(Sys.getenv("NSLOTS"))/4);
if(is.na(n.cores))n.cores <- 1;



qn.lower <- (l.trim/100-l.thres)/by.thres+1;
qn.upper <- (u.trim/100-l.thres)/by.thres+1;

n.delta <- NCOL(model.matrix(rho.formula, mysample));
rho.name <- paste("delta.",colnames(model.matrix(rho.formula, mysample)),sep="");


current.na.action <- options('na.action');
options(na.action='na.pass');
Xwc.all <- model.matrix(outcome.formula, mysample);
Zwc.all <- model.matrix(selection.formula, mysample);
X.name  <- colnames(Xwc.all);
n.par <- length(X.name);

options(current.na.action);


theta.counter1 <- elist.counter1$beta[,c(X.name,rho.name)];
theta.counter2 <- elist.counter2$beta[,c(X.name,rho.name)]; # n.thres-by-(K+3) matrix

pi.counter1 <- elist.counter1$pi[1,];
pi.counter2 <- elist.counter2$pi[1,];


F1   <- sort(apply(theta.counter1[,X.name], 1, Fy, X.counter=Xwc.all[mysample[,counter.id]==1,]));
F2   <- sort(apply(theta.counter2[,X.name], 1, Fy, X.counter=Xwc.all[mysample[,counter.id]==2,]));
F2_c <- sort(apply(theta.counter1[,X.name], 1, Fy, X.counter=Xwc.all[mysample[,counter.id]==2,]));
F1_c <- sort(apply(theta.counter2[,X.name], 1, Fy, X.counter=Xwc.all[mysample[,counter.id]==1,]));

QF1   <- approxfun(F1, y=ys, rule=2);
QF2_c <- approxfun(F2_c, y=ys, rule=2);
QF2   <- approxfun(F2, y=ys, rule=2);
pct.F.theta <- function(x){
  pct <- (QF2_c(x) - QF2(x))/(QF1(x) - QF2(x))*100;
#   pct[pct>100] <- 100;
#   pct[pct< -100] <- -100;
  return(pct)};
pct.DF.theta <- (F2_c - F2)/(F1 - F2)*100;
# pct.DF.theta[pct.DF.theta>100] <- 100;
# pct.DF.theta[pct.DF.theta< -100] <- -100;

e.avg.F1   <- ys[qn.upper] - integrate(approxfun(F1, x=ys, rule=2:1), ys[qn.lower], ys[qn.upper])$value;
e.avg.F2   <- ys[qn.upper] - integrate(approxfun(F2, x=ys, rule=2:1), ys[qn.lower], ys[qn.upper])$value;
e.avg.F2_c <- ys[qn.upper] - integrate(approxfun(F2_c, x=ys, rule=2:1), ys[qn.lower], ys[qn.upper])$value;
e.avg.F1_c <- ys[qn.upper] - integrate(approxfun(F1_c, x=ys, rule=2:1), ys[qn.lower], ys[qn.upper])$value;

P1 <- mean(pnorm(Zwc.all[mysample[,counter.id]==1,] %*% pi.counter1));
P2 <- mean(pnorm(Zwc.all[mysample[,counter.id]==2,] %*% pi.counter2));
P2_c <- mean(pnorm(Zwc.all[mysample[,counter.id]==2,] %*% pi.counter1));
P1_c <- mean(pnorm(Zwc.all[mysample[,counter.id]==1,] %*% pi.counter2));

beta.H.c1 <- cbind(ys-beta.Heckman.counter1[1], -t(replicate(n.thres,beta.Heckman.counter1[-1])))/sigma.Heckman.counter1;
beta.H.c2 <- cbind(ys-beta.Heckman.counter2[1], -t(replicate(n.thres,beta.Heckman.counter2[-1])))/sigma.Heckman.counter2;
F.Heck1   <- sort(apply(beta.H.c1, 1, Fy, X.counter=Xwc.all[mysample[,counter.id]==1,]));
F.Heck2   <- sort(apply(beta.H.c2, 1, Fy, X.counter=Xwc.all[mysample[,counter.id]==2,]));
F.Heck2_c <- sort(apply(beta.H.c1, 1, Fy, X.counter=Xwc.all[mysample[,counter.id]==2,]));
F.Heck1_c <- sort(apply(beta.H.c2, 1, Fy, X.counter=Xwc.all[mysample[,counter.id]==1,]));

QF.Heck1   <- approxfun(F.Heck1, y=ys, rule=2);
QF.Heck2_c <- approxfun(F.Heck2_c, y=ys, rule=2);
QF.Heck2   <- approxfun(F.Heck2, y=ys, rule=2);
pct.Heck.F.theta <- function(x){
  pct <- (QF.Heck2_c(x) - QF.Heck2(x))/(QF.Heck1(x) - QF.Heck2(x))*100;
#   pct[pct>100] <- 100;
#   pct[pct< -100] <- -100;
  return(pct)};
pct.Heck.DF.theta <- (F.Heck2_c - F.Heck2)/(F.Heck1 - F.Heck2)*100;
# pct.Heck.DF.theta[pct.Heck.DF.theta>100] <- 100;
# pct.Heck.DF.theta[pct.Heck.DF.theta< -100] <- -100;

rm(Xwc.all,Zwc.all);
##################################################
# Estimate the influence function and std errors #
##################################################

# 1. Counterfactual level 1 #
phi.para.1l <- mclapply(as.list(ys), phi.theta, mc.cores=n.cores,
                        parameter.ys=theta.counter1, ys.ind=ys, pi=pi.counter1, data.local=mysample, 
                        outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula,
                        counter.id=counter.id, population=1); # list of n.thres, each N-by-(K+1) vector 

phi.F_hat.1l <- mclapply(as.list(ys), phi.F, mc.cores=n.cores,
                         parameter.ys=theta.counter1, ys.ind=ys, F_hat.ys=F1, phi.para.ys=phi.para.1l, data.local=mysample, 
                         outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula,
                         counter.id=counter.id, population=1);

phi.F_hat.2_cl <- mclapply(as.list(ys), phi.F, mc.cores=n.cores,
                           parameter.ys=theta.counter1, ys.ind=ys, F_hat.ys=F2_c, phi.para.ys=phi.para.1l, data.local=mysample, 
                           outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula,
                           counter.id=counter.id, population=2);
# 
# # N*(n.par+n.delta)*n.thres, e.g. (i,k,y) means (observation,parameter,threshold) #
phi.para.1    <- do.call(abind, list(phi.para.1l,along=3));
phi.F_hat.1   <- do.call(cbind, phi.F_hat.1l); # N*n.thres
phi.F_hat.2_c <- do.call(cbind, phi.F_hat.2_cl);

phi.avg.1   <- sapply(apply(phi.F_hat.1, 1, approxfun, x=ys, rule=2:1),
                      function(myfunc){-integrate(myfunc,ys[qn.lower],ys[qn.upper],subdivisions=2000,stop.on.error=F)$value}); # N*1;
phi.avg.2_c <- sapply(apply(phi.F_hat.2_c, 1, approxfun, x=ys, rule=2:1),
                      function(myfunc){-integrate(myfunc,ys[qn.lower],ys[qn.upper],subdivisions=2000,stop.on.error=F)$value});


var.counter1 <- rbind(colMeans(phi.para.1^2), 
                      colMeans(phi.F_hat.1^2), 
                      colMeans(phi.F_hat.2_c^2))/NROW(phi.F_hat.1); # (K+3)*n.thres
rownames(var.counter1)[-(1:(n.par+n.delta))] <- c("F","F_counter");

var.avg.counter1 <- colMeans(cbind(phi.avg.1,phi.avg.2_c)^2)/NROW(phi.avg.1);


phi.counter1 <- abind(phi.para.1, phi.F_hat.1, phi.F_hat.2_c, along=2); # N*(n.par+n.delta+2)*n.thres array;

rm(phi.para.1l, phi.F_hat.1l, phi.F_hat.2_cl, phi.para.1);

###################################
# Bootstrap theta, F and Averages #
###################################
print(mem_used());
boot.counter1.l <- mclapply(1:n.boot, boot.func, mc.cores=n.cores, phi=phi.counter1, variance=var.counter1, 
                            var.boot=0);

rm(phi.counter1);

boot.avg.counter1.l <- mclapply(1:n.boot, boot.func, mc.cores=n.cores, phi=cbind(phi.avg.1,phi.avg.2_c), 
                                variance=var.avg.counter1, var.boot=0);

rm(phi.avg.1, phi.avg.2_c);





# 2. Counterfactual level 2 #

phi.para.2l <- mclapply(as.list(ys), phi.theta, mc.cores=n.cores,
                        parameter.ys=theta.counter2, ys.ind=ys, pi=pi.counter2, data.local=mysample, 
                        outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula,
                        counter.id=counter.id, population=2);

phi.F_hat.2l <- mclapply(as.list(ys), phi.F, mc.cores=n.cores,
                         parameter.ys=theta.counter2, ys.ind=ys, F_hat.ys=F2, phi.para.ys=phi.para.2l, data.local=mysample, 
                         outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula,
                         counter.id=counter.id, population=2);

phi.F_hat.1_cl <- mclapply(as.list(ys), phi.F, mc.cores=n.cores,
                           parameter.ys=theta.counter2, ys.ind=ys, F_hat.ys=F1_c, phi.para.ys=phi.para.2l, data.local=mysample, 
                           outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula,
                           counter.id=counter.id, population=1);

# N*(n.par+n.delta)*n.thres, e.g. (i,k,y) means (observation,parameter,threshold) #
phi.para.2    <- do.call(abind, list(phi.para.2l,along=3)); 
phi.F_hat.2   <- do.call(cbind, phi.F_hat.2l); # N*n.thres
phi.F_hat.1_c <- do.call(cbind, phi.F_hat.1_cl);

phi.avg.2   <- sapply(apply(phi.F_hat.2, 1, approxfun, x=ys, rule=2:1),
                      function(myfunc){-integrate(myfunc,ys[qn.lower],ys[qn.upper],subdivisions=2000,stop.on.error=F)$value}); # N*1;
phi.avg.1_c <- sapply(apply(phi.F_hat.1_c, 1, approxfun, x=ys, rule=2:1),
                      function(myfunc){-integrate(myfunc,ys[qn.lower],ys[qn.upper],subdivisions=2000,stop.on.error=F)$value});



var.counter2 <- rbind(colMeans(phi.para.2^2), 
                      colMeans(phi.F_hat.2^2), 
                      colMeans(phi.F_hat.1_c^2))/NROW(phi.F_hat.2);
rownames(var.counter2)[-(1:(n.par+n.delta))] <- c("F","F_counter");

var.avg.counter2 <- colMeans(cbind(phi.avg.2,phi.avg.1_c)^2)/NROW(phi.avg.2);


boot.pct.F.theta.l <- mclapply(1:n.boot, boot.func.pct, mc.cores=n.cores, F.tl=F2_c, F.tr=F2, F.bl=F1, F.br=F2, 
                               phi.tl=phi.F_hat.2_c, phi.tr=phi.F_hat.2, phi.bl=phi.F_hat.1, phi.br=phi.F_hat.2, ys=ys,
                               taus=seq(l.thres,u.thres,by.thres));
boot.pct.DF.theta.l <- mclapply(1:n.boot, boot.func.pct.DF, mc.cores=n.cores, F.tl=F2_c, F.tr=F2, F.bl=F1, F.br=F2, 
                                phi.tl=phi.F_hat.2_c, phi.tr=phi.F_hat.2, phi.bl=phi.F_hat.1, phi.br=phi.F_hat.2);
boot.QF.l <- mclapply(1:n.boot, boot.func.QF, mc.cores=n.cores, F.hat=cbind(F1,F2_c,F2,F1_c), 
                      phi.F=cbind(phi.F_hat.1, phi.F_hat.2_c, phi.F_hat.2, phi.F_hat.1_c), ys=ys, taus=seq(l.thres,u.thres,by.thres));


phi.counter2 <- abind(phi.para.2, phi.F_hat.2, phi.F_hat.1_c, along=2); # N*(n.par+n.delta+2)*n.thres array;

rm(phi.para.2l, phi.F_hat.2l, phi.F_hat.1_cl, phi.para.2,
   phi.F_hat.2, phi.F_hat.1_c, phi.F_hat.1, phi.F_hat.2_c);



###################################
# Bootstrap theta, F and Averages #
###################################
print(mem_used());
boot.counter2.l <- mclapply(1:n.boot, boot.func, mc.cores=n.cores, phi=phi.counter2, variance=var.counter2, 
                            var.boot=0);

rm(phi.counter2);

boot.avg.counter2.l <- mclapply(1:n.boot, boot.func, mc.cores=n.cores, phi=cbind(phi.avg.2,phi.avg.1_c), 
                                variance=var.avg.counter2, var.boot=0);

rm(phi.avg.2, phi.avg.1_c);


####################################
# Bootstrap for Selection Equation #
####################################
phi.P.1 <- phi.P(pi=pi.counter1, P=P1, data.local=mysample, outcome.formula=outcome.cormula, 
                 selection.formula=selection.formula, rho.formula=rho.formula, counter.id=counter.id, population=1);
phi.P.2_c <- phi.P(pi=pi.counter1, P=P2_c, data.local=mysample, outcome.formula=outcome.cormula, 
                   selection.formula=selection.formula, rho.formula=rho.formula, counter.id=counter.id, population=2);
phi.P.2 <- phi.P(pi=pi.counter2, P=P2, data.local=mysample, outcome.formula=outcome.cormula, 
                 selection.formula=selection.formula, rho.formula=rho.formula, counter.id=counter.id, population=2);
phi.P.1_c <- phi.P(pi=pi.counter2, P=P1_c, data.local=mysample, outcome.formula=outcome.cormula, 
                   selection.formula=selection.formula, rho.formula=rho.formula, counter.id=counter.id, population=1);

var.P <- colMeans(cbind(phi.P.1,phi.P.2_c,phi.P.2,phi.P.1_c)^2)/NROW(phi.P.1);

boot.P.l <- mclapply(1:n.boot, boot.func, mc.cores=n.cores, phi=cbind(phi.P.1,phi.P.2_c,phi.P.2,phi.P.1_c), 
                     variance=var.P, var.boot=0);

############################
# Critical Value and CB/CI #
############################
boot.counter1 <- do.call(abind,list(boot.counter1.l,along=3)); # (n.thres+n.thres)*(n.par+n.delta+2)*B matrix -- (thres,para,b)
boot.counter2 <- do.call(abind,list(boot.counter2.l,along=3)); 

delta.counter1 <- boot.counter1[1:n.thres, ,]; # n.thres*(n.par+n.delta+2)*B
delta.counter2 <- boot.counter2[1:n.thres, ,];

t.counter1 <- boot.counter1[(n.thres+1):(2*n.thres), ,]; # n.thres*(n.par+n.delta+2)*B
t.counter2 <- boot.counter2[(n.thres+1):(2*n.thres), ,];


boot.avg.counter1 <- do.call(abind,list(boot.avg.counter1.l,along=3));
boot.avg.counter2 <- do.call(abind,list(boot.avg.counter2.l,along=3));

t.avg.counter1 <- boot.avg.counter1[2,,];
t.avg.counter2 <- boot.avg.counter2[2,,];

boot.P <- do.call(abind,list(boot.P.l,along=3));
t.P <- boot.P[2,,];

boot.pct.F.theta <- do.call(cbind,boot.pct.F.theta.l);
delta.pct.F.theta <- boot.pct.F.theta - replicate(n.boot,pct.F.theta(seq(l.thres,u.thres,by.thres)));
se.pct.F.theta <- apply(delta.pct.F.theta,1,sd,na.rm=T);
t.pct.F.theta <- abs(delta.pct.F.theta)/se.pct.F.theta;

t.F.joint <- apply(abind(array(t.counter1[,n.par+n.delta+1,],c(n.thres,1,n.boot)),
                         t.counter2[,n.par+n.delta+1,],
                         t.counter1[,n.par+n.delta+2,],along=2), c(1,3), max, na.rm=T);


boot.QF <- do.call(abind, list(boot.QF.l, along=3));
delta.QF <- boot.QF[,1:3,] - replicate(n.boot, cbind(QF1(seq(l.thres,u.thres,by.thres)),
                                                     QF2_c(seq(l.thres,u.thres,by.thres)),
                                                     QF2(seq(l.thres,u.thres,by.thres))));
se.QF <- apply(delta.QF, c(1,2), sd, na.rm=T);
t.QF.l <- lapply(1:n.boot, function(b){abs(delta.QF[,,b])/se.QF});
t.QF <- do.call(abind, list(t.QF.l, along=3)); # n.thres*4*n.boot;
t.QF.joint <- apply(t.QF, c(1,3), max, na.rm=T); # n.thres*n.boot;


boot.pct.DF.theta  <- do.call(cbind, boot.pct.DF.theta.l);
delta.pct.DF.theta <- boot.pct.DF.theta - replicate(n.boot,pct.DF.theta);
se.pct.DF.theta    <- apply(delta.pct.DF.theta,1,sd,na.rm=T);
t.pct.DF.theta     <- abs(delta.pct.DF.theta)/se.pct.DF.theta;



crt.counter1 <- apply(apply(t.counter1[qn.lower:qn.upper,,],c(2:3),max,na.rm=T),1,quantile,1-alpha,na.rm=T); # (n.par+n.delta+2)*1
crt.counter2 <- apply(apply(t.counter2[qn.lower:qn.upper,,],c(2:3),max,na.rm=T),1,quantile,1-alpha,na.rm=T);

crt.avg.counter1 <- apply(t.avg.counter1, 1, quantile, 1-alpha, na.rm=T);
crt.avg.counter2 <- apply(t.avg.counter2, 1, quantile, 1-alpha, na.rm=T);

crt.P <- apply(t.P, 1, quantile, 1-alpha, na.rm=T);

crt.pct.F.theta <- quantile(apply(t.pct.F.theta[qn.lower:qn.upper,], 2, max, na.rm=T), 1-alpha, na.rm=T);
crt.F.joint <- quantile(apply(t.F.joint[qn.lower:qn.upper,],2,max,na.rm=T), 1-alpha, na.rm=T);

crt.avg.joint <- quantile(apply(rbind(t.avg.counter1,t.avg.counter2[1,]),2,max,na.rm=T), 1-alpha, na.rm=T);
crt.P.joint <- quantile(apply(t.P[1:3,],2,max,na.rm=T), 1-alpha, na.rm=T); 

crt.QF.joint <- quantile(apply(t.QF.joint[qn.lower:qn.upper,],2,max,na.rm=T), 1-alpha, na.rm=T);
crt.pct.DF.theta <- quantile(apply(t.pct.DF.theta[qn.lower:qn.upper,],2,max,na.rm=T), 1-alpha, na.rm=T);


ubound.counter1 <- cbind(theta.counter1,F1,F2_c) + t(crt.counter1*sqrt(var.counter1));
lbound.counter1 <- cbind(theta.counter1,F1,F2_c) - t(crt.counter1*sqrt(var.counter1));

ubound.counter2 <- cbind(theta.counter2,F2,F1_c) + t(crt.counter2*sqrt(var.counter2));
lbound.counter2 <- cbind(theta.counter2,F2,F1_c) - t(crt.counter2*sqrt(var.counter2));

ubound.F1 <- sort(ubound.counter1[,n.par+n.delta+1]);
lbound.F1 <- sort(lbound.counter1[,n.par+n.delta+1]);

ubound.F2 <- sort(ubound.counter2[,n.par+n.delta+1]);
lbound.F2 <- sort(lbound.counter2[,n.par+n.delta+1]);

ubound.F2_c <- sort(ubound.counter1[,n.par+n.delta+2]);
lbound.F2_c <- sort(lbound.counter1[,n.par+n.delta+2]);

ubound.F1_c <- sort(ubound.counter2[,n.par+n.delta+2]);
lbound.F1_c <- sort(lbound.counter2[,n.par+n.delta+2]);


ubound.F1[ubound.F1 > 1] <- 1;
lbound.F1[lbound.F1 < 0] <- 0;

ubound.F2[ubound.F2 > 1] <- 1;
lbound.F2[lbound.F2 < 0] <- 0;

ubound.F1_c[ubound.F1_c > 1] <- 1;
lbound.F1_c[lbound.F1_c < 0] <- 0;

ubound.F2_c[ubound.F2_c > 1] <- 1;
lbound.F2_c[lbound.F2_c < 0] <- 0;


ubound.avg <- c(e.avg.F1, e.avg.F2_c, e.avg.F2, e.avg.F1_c) + c(crt.avg.counter1*sqrt(var.avg.counter1),
                                                                crt.avg.counter2*sqrt(var.avg.counter2));
lbound.avg <- c(e.avg.F1, e.avg.F2_c, e.avg.F2, e.avg.F1_c) - c(crt.avg.counter1*sqrt(var.avg.counter1),
                                                                crt.avg.counter2*sqrt(var.avg.counter2));
ubound.avg.joint <- c(e.avg.F1, e.avg.F2_c, e.avg.F2) + crt.avg.joint*sqrt(c(var.avg.counter1,var.avg.counter2[1]));
lbound.avg.joint <- c(e.avg.F1, e.avg.F2_c, e.avg.F2) - crt.avg.joint*sqrt(c(var.avg.counter1,var.avg.counter2[1]));

ubound.P <- c(P1, P2_c, P2, P1_c) + crt.P*sqrt(var.P);
lbound.P <- c(P1, P2_c, P2, P1_c) - crt.P*sqrt(var.P);
ubound.P.joint <- c(P1, P2_c, P2) + crt.P.joint*sqrt(var.P[1:3]);
lbound.P.joint <- c(P1, P2_c, P2) - crt.P.joint*sqrt(var.P[1:3]);

ubound.pct.F.theta <- function(x){
  u.pct <- pct.F.theta(x) + crt.pct.F.theta*approxfun(se.pct.F.theta,x=seq(l.thres,u.thres,by.thres),rule=2)(x);
#   u.pct[u.pct>100] <- 100;
#   u.pct[u.pct< -100] <- -100;
  return(u.pct)};
lbound.pct.F.theta <- function(x){
  l.pct <- pct.F.theta(x) - crt.pct.F.theta*approxfun(se.pct.F.theta,x=seq(l.thres,u.thres,by.thres),rule=2)(x);
#   l.pct[l.pct>100] <- 100;
#   l.pct[l.pct< -100] <- -100;
  return(l.pct)};

ubound.F.joint <- apply(cbind(F1,F2_c,F2)+t(crt.F.joint*
                                              sqrt(rbind(var.counter1[c("F","F_counter"),],var.counter2["F",]))),2,sort);
lbound.F.joint <- apply(cbind(F1,F2_c,F2)-t(crt.F.joint*
                                              sqrt(rbind(var.counter1[c("F","F_counter"),],var.counter2["F",]))),2,sort);
ubound.F.joint[ubound.F.joint>1] <- 1;
ubound.F.joint[ubound.F.joint<0] <- 0;
lbound.F.joint[lbound.F.joint>1] <- 1;
lbound.F.joint[lbound.F.joint<0] <- 0;
ubound.QF.joint <- apply(ubound.F.joint, 2, approxfun, y=ys, rule=2);
lbound.QF.joint <- apply(lbound.F.joint, 2, approxfun, y=ys, rule=2);

ubound.QF <- cbind(QF1(seq(l.thres,u.thres,by.thres)),
                   QF2_c(seq(l.thres,u.thres,by.thres)),
                   QF2(seq(l.thres,u.thres,by.thres))) + crt.QF.joint*se.QF;
lbound.QF <- cbind(QF1(seq(l.thres,u.thres,by.thres)),
                   QF2_c(seq(l.thres,u.thres,by.thres)),
                   QF2(seq(l.thres,u.thres,by.thres))) - crt.QF.joint*se.QF;

ubound.pct.DF.theta <- pct.DF.theta + crt.pct.DF.theta*se.pct.DF.theta;
lbound.pct.DF.theta <- pct.DF.theta - crt.pct.DF.theta*se.pct.DF.theta;

# ubound.pct.DF.theta[ubound.pct.DF.theta>100] <- 100;
# ubound.pct.DF.theta[ubound.pct.DF.theta< -100] <- -100;
# lbound.pct.DF.theta[lbound.pct.DF.theta>100] <- 100;
# lbound.pct.DF.theta[lbound.pct.DF.theta< -100] <- 100;

################
# Save Results #
################
save.image(paste("Results",str.sample,str.rho,str.trim,str.cb,str.robust,str.boot,".RData",sep=""));

# load(paste("Results",str.sample,str.rho,str.trim,".RData",sep=""));
# source("DRS_Source.R");
# # Calculate the namber of cores
# n.cores <- ceiling(as.numeric(Sys.getenv("NSLOTS"))/4);
# if(is.na(n.cores))n.cores <- 1;

#######################################################################################################################
#################################################
# Empirical Distribution of Observed Population #
#################################################
F1.obs <- sapply(ys,function(y){mean((mysample[mysample[,counter.id]==1, Y.name]<=y),na.rm=T)});
F2.obs <- sapply(ys,function(y){mean((mysample[mysample[,counter.id]==2, Y.name]<=y),na.rm=T)});

##############################################
# Fitted Distribution of Observed Population #
##############################################
F2.condl <- mclapply(1:n.thres, Fy.cond, mc.cores=n.cores,
                     pi=pi.counter2, beta=theta.counter2[,X.name], athrho=theta.counter2[,rho.name],
                     sample.counter=mysample[mysample[,counter.id]==2,],
                     outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula);
F2.cond <- sort(do.call(c,F2.condl));
F2.cond.c1rhol <- mclapply(1:n.thres,Fy.cond, mc.cores=n.cores,
                           pi=pi.counter2, beta=theta.counter2[,X.name],athrho=theta.counter1[,rho.name],
                           sample.counter=mysample[mysample[,counter.id]==2,],
                           outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula);
F2.cond.c1rho <- sort(do.call(c,F2.cond.c1rhol));
F2.cond.c1rhopil <- mclapply(1:n.thres,Fy.cond, mc.cores=n.cores,
                             pi=pi.counter1, beta=theta.counter2[,X.name],athrho=theta.counter1[,rho.name],
                             sample.counter=mysample[mysample[,counter.id]==2,],
                             outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula);
F2.cond.c1rhopi <- sort(do.call(c,F2.cond.c1rhopil));
F2.cond.c1rhopibetal <- mclapply(1:n.thres,Fy.cond, mc.cores=n.cores,
                                 pi=pi.counter1, beta=theta.counter1[,X.name],athrho=theta.counter1[,rho.name],
                                 sample.counter=mysample[mysample[,counter.id]==2,],
                                 outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula);
F2.cond.c1rhopibeta <- sort(do.call(c,F2.cond.c1rhopibetal));
F1.condl <- mclapply(1:n.thres,Fy.cond, mc.cores=n.cores,
                     pi=pi.counter1,beta=theta.counter1[,X.name],athrho=theta.counter1[,rho.name],
                     sample.counter=mysample[mysample[,counter.id]==1,],
                     outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula);
F1.cond <- sort(do.call(c,F1.condl));

QF2.cond             <- approxfun(F2.cond, y=ys, rule=2);
QF2.cond.c1rho       <- approxfun(F2.cond.c1rho, y=ys, rule=2);
QF2.cond.c1rhopi     <- approxfun(F2.cond.c1rhopi, y=ys, rule=2);
QF2.cond.c1rhopibeta <- approxfun(F2.cond.c1rhopibeta, y=ys, rule=2);
QF1.cond             <- approxfun(F1.cond, y=ys, rule=2);
pct.Fcond.rho  <- function(x){(QF2.cond.c1rho(x) - QF2.cond(x))/(QF1.cond(x) - QF2.cond(x))*100};
pct.Fcond.pi   <- function(x){(QF2.cond.c1rhopi(x) - QF2.cond.c1rho(x))/(QF1.cond(x) - QF2.cond(x))*100};
pct.Fcond.beta <- function(x){(QF2.cond.c1rhopibeta(x) - QF2.cond.c1rhopi(x))/(QF1.cond(x) - QF2.cond(x))*100};
pct.Fcond.Fz   <- function(x){(QF1.cond(x) - QF2.cond.c1rhopibeta(x))/(QF1.cond(x) - QF2.cond(x))*100};
pct.Fcond.selection <- function(x){(QF2.cond.c1rhopi(x) - QF2.cond(x))/(QF1.cond(x) - QF2.cond(x))*100};

pct.DFcond.rho  <- (F2.cond.c1rho - F2.cond)/(F1.cond - F2.cond)*100;
pct.DFcond.pi   <- (F2.cond.c1rhopi - F2.cond.c1rho)/(F1.cond - F2.cond)*100;
pct.DFcond.beta <- (F2.cond.c1rhopibeta - F2.cond.c1rhopi)/(F1.cond - F2.cond)*100;
pct.DFcond.Fz   <- (F1.cond - F2.cond.c1rhopibeta)/(F1.cond - F2.cond)*100;
pct.DFcond.selection <- (F2.cond.c1rhopi - F2.cond)/(F1.cond - F2.cond)*100;

e.avg.F2.cond             <- ys[qn.upper] - integrate(approxfun(F2.cond, x=ys, rule=2:1), ys[qn.lower], ys[qn.upper])$value;
e.avg.F2.cond.c1rho       <- ys[qn.upper] - integrate(approxfun(F2.cond.c1rho, x=ys, rule=2:1), ys[qn.lower], ys[qn.upper])$value;
e.avg.F2.cond.c1rhopi     <- ys[qn.upper] - integrate(approxfun(F2.cond.c1rhopi, x=ys, rule=2:1), ys[qn.lower], ys[qn.upper])$value;
e.avg.F2.cond.c1rhopibeta <- ys[qn.upper] - integrate(approxfun(F2.cond.c1rhopibeta, x=ys, rule=2:1), ys[qn.lower], ys[qn.upper])$value;
e.avg.F1.cond             <- ys[qn.upper] - integrate(approxfun(F1.cond, x=ys, rule=2:1), ys[qn.lower], ys[qn.upper])$value;

F2.Heck.condl <- mclapply(1:n.thres, Fy.cond, mc.cores=n.cores,
                          pi=pi.Heckman.counter2, beta=beta.H.c2, athrho=-rho.Heckman.counter2,
                          sample.counter=mysample[mysample[,counter.id]==2,],
                          outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=work~1);
F2.Heck.cond <- sort(do.call(c,F2.Heck.condl));
F2.Heck.cond.c1rhol <- mclapply(1:n.thres,Fy.cond, mc.cores=n.cores,
                                pi=pi.Heckman.counter2, beta=beta.H.c2, athrho=-rho.Heckman.counter1,
                                sample.counter=mysample[mysample[,counter.id]==2,],
                                outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=work~1);
F2.Heck.cond.c1rho <- sort(do.call(c,F2.Heck.cond.c1rhol));
F2.Heck.cond.c1rhopil <- mclapply(1:n.thres,Fy.cond, mc.cores=n.cores,
                                  pi=pi.Heckman.counter1, beta=beta.H.c2, athrho=-rho.Heckman.counter1,
                                  sample.counter=mysample[mysample[,counter.id]==2,],
                                  outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=work~1);
F2.Heck.cond.c1rhopi <- sort(do.call(c,F2.Heck.cond.c1rhopil));
F2.Heck.cond.c1rhopibetal <- mclapply(1:n.thres,Fy.cond, mc.cores=n.cores,
                                      pi=pi.Heckman.counter1, beta=beta.H.c1, athrho=-rho.Heckman.counter1,
                                      sample.counter=mysample[mysample[,counter.id]==2,],
                                      outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=work~1);
F2.Heck.cond.c1rhopibeta <- sort(do.call(c,F2.Heck.cond.c1rhopibetal));
F1.Heck.condl <- mclapply(1:n.thres,Fy.cond, mc.cores=n.cores,
                          pi=pi.Heckman.counter1, beta=beta.H.c1, athrho=-rho.Heckman.counter1,
                          sample.counter=mysample[mysample[,counter.id]==1,],
                          outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=work~1);
F1.Heck.cond <- sort(do.call(c,F1.Heck.condl));


QF2.Heck.cond             <- approxfun(F2.Heck.cond, y=ys, rule=2);
QF2.Heck.cond.c1rho       <- approxfun(F2.Heck.cond.c1rho, y=ys, rule=2);
QF2.Heck.cond.c1rhopi     <- approxfun(F2.Heck.cond.c1rhopi, y=ys, rule=2);
QF2.Heck.cond.c1rhopibeta <- approxfun(F2.Heck.cond.c1rhopibeta, y=ys, rule=2);
QF1.Heck.cond             <- approxfun(F1.Heck.cond, y=ys, rule=2);
pct.Heck.Fcond.rho  <- function(x){(QF2.Heck.cond.c1rho(x) - QF2.Heck.cond(x))/(QF1.Heck.cond(x) - QF2.Heck.cond(x))*100};
pct.Heck.Fcond.pi   <- function(x){(QF2.Heck.cond.c1rhopi(x) - QF2.Heck.cond.c1rho(x))/(QF1.Heck.cond(x) - QF2.Heck.cond(x))*100};
pct.Heck.Fcond.beta <- function(x){(QF2.Heck.cond.c1rhopibeta(x) - QF2.Heck.cond.c1rhopi(x))/(QF1.Heck.cond(x) - QF2.Heck.cond(x))*100};
pct.Heck.Fcond.Fz   <- function(x){(QF1.Heck.cond(x) - QF2.Heck.cond.c1rhopibeta(x))/(QF1.Heck.cond(x) - QF2.Heck.cond(x))*100};
pct.Heck.Fcond.selection <- function(x){(QF2.Heck.cond.c1rhopi(x) - QF2.Heck.cond(x))/(QF1.Heck.cond(x) - QF2.Heck.cond(x))*100};

pct.Heck.DFcond.rho  <- (F2.Heck.cond.c1rho - F2.Heck.cond)/(F1.Heck.cond - F2.Heck.cond)*100;
pct.Heck.DFcond.pi   <- (F2.Heck.cond.c1rhopi - F2.Heck.cond.c1rho)/(F1.Heck.cond - F2.Heck.cond)*100;
pct.Heck.DFcond.beta <- (F2.Heck.cond.c1rhopibeta - F2.Heck.cond.c1rhopi)/(F1.Heck.cond - F2.Heck.cond)*100;
pct.Heck.DFcond.Fz   <- (F1.Heck.cond - F2.Heck.cond.c1rhopibeta)/(F1.Heck.cond - F2.Heck.cond)*100;
pct.Heck.DFcond.selection <- (F2.Heck.cond.c1rhopi - F2.Heck.cond)/(F1.Heck.cond - F2.Heck.cond)*100;

##########################
## Setting rho = 0 plot ##
##########################
F1.noselectionl <- mclapply(1:n.thres,Fy.cond, mc.cores=n.cores,
                            pi=pi.counter1,beta=theta.counter1[,X.name],athrho=0,
                            sample.counter=mysample[mysample[,counter.id]==1,],
                            outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula);
F1.noselection <- sort(do.call(c,F1.noselectionl));

F2.noselectionl <- mclapply(1:n.thres, Fy.cond, mc.cores=n.cores,
                            pi=pi.counter2, beta=theta.counter2[,X.name], athrho=0,
                            sample.counter=mysample[mysample[,counter.id]==2,],
                            outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula);
F2.noselection <- sort(do.call(c,F2.noselectionl));

##################################################
# Estimate the influence function and std errors #
##################################################
phi.pi.counter1 <- phi.pi(pi=pi.counter1, data.local=mysample, selection.formula, counter.id, population=1);
phi.pi.counter2 <- phi.pi(pi=pi.counter2, data.local=mysample, selection.formula, counter.id, population=2);
phi.para.counter1l <- mclapply(as.list(ys),phi.theta, mc.cores=n.cores,
                               parameter.ys=theta.counter1, ys.ind=ys, pi=pi.counter1, data.local=mysample, 
                               outcome.formula = outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula,
                               counter.id=counter.id, population=1); # list of n.thres, each N-by-(K+1) vector 
phi.para.counter2l <- mclapply(as.list(ys),phi.theta, mc.cores=n.cores,
                               parameter.ys=theta.counter2, ys.ind=ys, pi=pi.counter2, data.local=mysample, 
                               outcome.formula = outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula,
                               counter.id=counter.id, population=2); # list of n.thres, each N-by-(K+1) vector 

print(mem_used());
theta.all <- abind(theta.counter1,theta.counter2,along=3);
pi.all <- cbind(pi.counter1,pi.counter2);
# phi.para.*** or phi.pi.*** denote the influence functions of para(beta) or pi;
# phi.***.beta, phi.***.rho, phi.***.pi denote the three levels of counterfactual elements;
# population.v=c(Fz,beta,rho,pi);
phi.F2.condl <- mclapply(as.list(ys),phi.F.cond, mc.cores=n.cores,
                         parameter.ys=theta.all, ys.ind=ys, pi=pi.all, F_hat.ys=F2.cond, 
                         phi.para.ys.beta=phi.para.counter2l, phi.pi.beta=phi.pi.counter2, 
                         phi.para.ys.rho=phi.para.counter2l, phi.pi.rho=phi.pi.counter2, 
                         phi.pi.pi=phi.pi.counter2, data.local=mysample, 
                         outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula,
                         counter.id=counter.id, population.v=c(2,2,2,2));

phi.F2.cond.c1rhol <- mclapply(as.list(ys),phi.F.cond, mc.cores=n.cores,
                               parameter.ys=theta.all, ys.ind=ys, pi=pi.all, F_hat.ys=F2.cond.c1rho, 
                               phi.para.ys.beta=phi.para.counter2l, phi.pi.beta=phi.pi.counter2, 
                               phi.para.ys.rho=phi.para.counter1l, phi.pi.rho=phi.pi.counter1, 
                               phi.pi.pi=phi.pi.counter2, data.local=mysample, 
                               outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula, 
                               counter.id=counter.id, population.v=c(2,2,1,2));

phi.F2.cond.c1rhopil <- mclapply(as.list(ys),phi.F.cond, mc.cores=n.cores,
                                 parameter.ys=theta.all, ys.ind=ys, pi=pi.all, F_hat.ys=F2.cond.c1rhopi, 
                                 phi.para.ys.beta=phi.para.counter2l, phi.pi.beta=phi.pi.counter2, 
                                 phi.para.ys.rho=phi.para.counter1l, phi.pi.rho=phi.pi.counter1, 
                                 phi.pi.pi=phi.pi.counter1, data.local=mysample, 
                                 outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula, 
                                 counter.id=counter.id, population.v=c(2,2,1,1));

phi.F2.cond.c1rhopibetal <- mclapply(as.list(ys),phi.F.cond, mc.cores=n.cores,
                                     parameter.ys=theta.all, ys.ind=ys, pi=pi.all, F_hat.ys=F2.cond.c1rhopibeta, 
                                     phi.para.ys.beta=phi.para.counter1l, phi.pi.beta=phi.pi.counter1, 
                                     phi.para.ys.rho=phi.para.counter1l, phi.pi.rho=phi.pi.counter1, 
                                     phi.pi.pi=phi.pi.counter1, data.local=mysample, 
                                     outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula,
                                     counter.id=counter.id, population.v=c(2,1,1,1));

phi.F1.condl <- mclapply(as.list(ys),phi.F.cond, mc.cores=n.cores,
                         parameter.ys=theta.all, ys.ind=ys, pi=pi.all, F_hat.ys=F1.cond,
                         phi.para.ys.beta=phi.para.counter1l, phi.pi.beta=phi.pi.counter1, 
                         phi.para.ys.rho=phi.para.counter1l, phi.pi.rho=phi.pi.counter1, 
                         phi.pi.pi=phi.pi.counter1, data.local=mysample, 
                         outcome.formula=outcome.formula, selection.formula=selection.formula, rho.formula=rho.formula,
                         counter.id=counter.id, population.v=c(1,1,1,1));

rm(phi.pi.counter1, phi.pi.counter2, phi.para.counter1l, phi.para.counter2l);

phi.F2.cond             <- do.call(cbind, phi.F2.condl); # N*n.thres
phi.F2.cond.c1rho       <- do.call(cbind, phi.F2.cond.c1rhol);
phi.F2.cond.c1rhopi     <- do.call(cbind, phi.F2.cond.c1rhopil);
phi.F2.cond.c1rhopibeta <- do.call(cbind, phi.F2.cond.c1rhopibetal);
phi.F1.cond             <- do.call(cbind, phi.F1.condl);

phi.avg.F2.cond             <- sapply(apply(phi.F2.cond, 1, approxfun, x=ys, rule=2:1),
                                      function(myfunc){-integrate(myfunc,ys[qn.lower],ys[qn.upper],subdivisions=2000,stop.on.error=F)$value}); # N*1;
phi.avg.F2.cond.c1rho       <- sapply(apply(phi.F2.cond.c1rho, 1, approxfun, x=ys, rule=2:1),
                                      function(myfunc){-integrate(myfunc,ys[qn.lower],ys[qn.upper],subdivisions=2000,stop.on.error=F)$value});
phi.avg.F2.cond.c1rhopi     <- sapply(apply(phi.F2.cond.c1rhopi, 1, approxfun, x=ys, rule=2:1),
                                      function(myfunc){-integrate(myfunc,ys[qn.lower],ys[qn.upper],subdivisions=2000,stop.on.error=F)$value});
phi.avg.F2.cond.c1rhopibeta <- sapply(apply(phi.F2.cond.c1rhopibeta, 1, approxfun, x=ys, rule=2:1),
                                      function(myfunc){-integrate(myfunc,ys[qn.lower],ys[qn.upper],subdivisions=2000,stop.on.error=F)$value});
phi.avg.F1.cond             <- sapply(apply(phi.F1.cond, 1, approxfun, x=ys, rule=2:1),
                                      function(myfunc){-integrate(myfunc,ys[qn.lower],ys[qn.upper],subdivisions=2000,stop.on.error=F)$value});


var.cond <- rbind(colMeans(phi.F2.cond^2), 
                  colMeans(phi.F2.cond.c1rho^2),
                  colMeans(phi.F2.cond.c1rhopi^2),
                  colMeans(phi.F2.cond.c1rhopibeta^2),
                  colMeans(phi.F1.cond^2))/NROW(phi.F2.cond); # 5*n.thres
rownames(var.cond) <- c("F2","F2|c1rho","F2|c1rhopi","F2|c1rhopibeta","F1");

var.avg.cond <- colMeans(cbind(phi.avg.F2.cond,
                               phi.avg.F2.cond.c1rho,
                               phi.avg.F2.cond.c1rhopi,
                               phi.avg.F2.cond.c1rhopibeta,
                               phi.avg.F1.cond)^2)/NROW(phi.avg.F2.cond);


phi.cond <- abind(array(phi.F2.cond,dim=c(NROW(phi.F2.cond),1,n.thres)), phi.F2.cond.c1rho, 
                  phi.F2.cond.c1rhopi, phi.F2.cond.c1rhopibeta, phi.F1.cond, along=2);
# N*5*n.thres array

boot.pct.Fcond.rho.l <- mclapply(1:n.boot, boot.func.pct, mc.cores=n.cores, F.tl=F2.cond.c1rho, F.tr=F2.cond, 
                                 F.bl=F1.cond, F.br=F2.cond, phi.tl=phi.F2.cond.c1rho, phi.tr=phi.F2.cond, 
                                 phi.bl=phi.F1.cond, phi.br=phi.F2.cond, ys=ys, taus=seq(l.thres,u.thres,by.thres));
boot.pct.Fcond.pi.l <- mclapply(1:n.boot, boot.func.pct, mc.cores=n.cores, F.tl=F2.cond.c1rhopi, F.tr=F2.cond.c1rho, 
                                F.bl=F1.cond, F.br=F2.cond, phi.tl=phi.F2.cond.c1rhopi, phi.tr=phi.F2.cond.c1rho, 
                                phi.bl=phi.F1.cond, phi.br=phi.F2.cond, ys=ys, taus=seq(l.thres,u.thres,by.thres));
boot.pct.Fcond.beta.l <- mclapply(1:n.boot, boot.func.pct, mc.cores=n.cores, F.tl=F2.cond.c1rhopibeta, F.tr=F2.cond.c1rhopi, 
                                  F.bl=F1.cond, F.br=F2.cond, phi.tl=phi.F2.cond.c1rhopibeta, phi.tr=phi.F2.cond.c1rhopi, 
                                  phi.bl=phi.F1.cond, phi.br=phi.F2.cond, ys=ys, taus=seq(l.thres,u.thres,by.thres));
boot.pct.Fcond.Fz.l <- mclapply(1:n.boot, boot.func.pct, mc.cores=n.cores, F.tl=F1.cond, F.tr=F2.cond.c1rhopibeta, 
                                F.bl=F1.cond, F.br=F2.cond, phi.tl=phi.F1.cond, phi.tr=phi.F2.cond.c1rhopibeta, 
                                phi.bl=phi.F1.cond, phi.br=phi.F2.cond, ys=ys, taus=seq(l.thres,u.thres,by.thres));
boot.pct.Fcond.selection.l <- mclapply(1:n.boot, boot.func.pct, mc.cores=n.cores, F.tl=F2.cond.c1rhopi, F.tr=F2.cond, 
                                  F.bl=F1.cond, F.br=F2.cond, phi.tl=phi.F2.cond.c1rhopi, phi.tr=phi.F2.cond, 
                                  phi.bl=phi.F1.cond, phi.br=phi.F2.cond, ys=ys, taus=seq(l.thres,u.thres,by.thres));

boot.pct.DFcond.rho.l <- mclapply(1:n.boot, boot.func.pct.DF, mc.cores=n.cores, F.tl=F2.cond.c1rho, F.tr=F2.cond, 
                                 F.bl=F1.cond, F.br=F2.cond, phi.tl=phi.F2.cond.c1rho, phi.tr=phi.F2.cond, 
                                 phi.bl=phi.F1.cond, phi.br=phi.F2.cond);
boot.pct.DFcond.pi.l <- mclapply(1:n.boot, boot.func.pct.DF, mc.cores=n.cores, F.tl=F2.cond.c1rhopi, F.tr=F2.cond.c1rho, 
                                F.bl=F1.cond, F.br=F2.cond, phi.tl=phi.F2.cond.c1rhopi, phi.tr=phi.F2.cond.c1rho, 
                                phi.bl=phi.F1.cond, phi.br=phi.F2.cond);
boot.pct.DFcond.beta.l <- mclapply(1:n.boot, boot.func.pct.DF, mc.cores=n.cores, F.tl=F2.cond.c1rhopibeta, F.tr=F2.cond.c1rhopi, 
                                  F.bl=F1.cond, F.br=F2.cond, phi.tl=phi.F2.cond.c1rhopibeta, phi.tr=phi.F2.cond.c1rhopi, 
                                  phi.bl=phi.F1.cond, phi.br=phi.F2.cond);
boot.pct.DFcond.Fz.l <- mclapply(1:n.boot, boot.func.pct.DF, mc.cores=n.cores, F.tl=F1.cond, F.tr=F2.cond.c1rhopibeta, 
                                F.bl=F1.cond, F.br=F2.cond, phi.tl=phi.F1.cond, phi.tr=phi.F2.cond.c1rhopibeta, 
                                phi.bl=phi.F1.cond, phi.br=phi.F2.cond);
boot.pct.DFcond.selection.l <- mclapply(1:n.boot, boot.func.pct.DF, mc.cores=n.cores, F.tl=F2.cond.c1rhopi, F.tr=F2.cond, 
                                       F.bl=F1.cond, F.br=F2.cond, phi.tl=phi.F2.cond.c1rhopi, phi.tr=phi.F2.cond, 
                                       phi.bl=phi.F1.cond, phi.br=phi.F2.cond);
boot.QFcond.l <- mclapply(1:n.boot, boot.func.QF, mc.cores=n.cores, 
                          F.hat=cbind(F2.cond, F2.cond.c1rho, F2.cond.c1rhopi, F2.cond.c1rhopibeta, F1.cond), 
                          phi.F=cbind(phi.F2.cond,phi.F2.cond.c1rho,phi.F2.cond.c1rhopi,phi.F2.cond.c1rhopibeta,phi.F1.cond),
                          ys=ys, taus=seq(l.thres,u.thres,by.thres));


rm(phi.F2.cond, phi.F2.cond.c1rho, phi.F2.cond.c1rhopi, phi.F2.cond.c1rhopibeta, phi.F1.cond,
   phi.F2.condl, phi.F2.cond.c1rhol, phi.F2.cond.c1rhopil, phi.F2.cond.c1rhopibetal, phi.F1.condl);



#########################
# Bootstrap theta and F #
#########################
print(mem_used());
boot.cond.l   <- mclapply(1:n.boot, boot.func, mc.cores=n.cores, phi = phi.cond, variance=var.cond, var.boot=0);
boot.cond <- do.call(abind,list(boot.cond.l,along=3)); # (n.thres+n.thres)*5*B matrix -- (thres,F,b)

rm(phi.cond, boot.cond.l);


boot.avg.cond.l <- mclapply(1:n.boot, boot.func, mc.cores=n.cores, phi=cbind(phi.avg.F2.cond,
                                                                             phi.avg.F2.cond.c1rho,
                                                                             phi.avg.F2.cond.c1rhopi,
                                                                             phi.avg.F2.cond.c1rhopibeta,
                                                                             phi.avg.F1.cond), variance=var.avg.cond, 
                            var.boot=0);
boot.avg.cond <- do.call(abind,list(boot.avg.cond.l,along=3));


boot.pct.Fcond.rho  <- do.call(cbind, boot.pct.Fcond.rho.l);
boot.pct.Fcond.pi   <- do.call(cbind, boot.pct.Fcond.pi.l);
boot.pct.Fcond.beta <- do.call(cbind, boot.pct.Fcond.beta.l);
boot.pct.Fcond.Fz   <- do.call(cbind, boot.pct.Fcond.Fz.l);
boot.pct.Fcond.selection   <- do.call(cbind, boot.pct.Fcond.selection.l);

boot.pct.DFcond.rho  <- do.call(cbind, boot.pct.DFcond.rho.l);
boot.pct.DFcond.pi   <- do.call(cbind, boot.pct.DFcond.pi.l);
boot.pct.DFcond.beta <- do.call(cbind, boot.pct.DFcond.beta.l);
boot.pct.DFcond.Fz   <- do.call(cbind, boot.pct.DFcond.Fz.l);
boot.pct.DFcond.selection   <- do.call(cbind, boot.pct.DFcond.selection.l);

boot.QFcond <- do.call(abind, list(boot.QFcond.l,along=3));

delta.cond <- boot.cond[1:n.thres, ,]; # n.thres*5*B

t.cond <- boot.cond[(n.thres+1):(2*n.thres), ,]; # n.thres*5*B
t.avg.cond <- boot.avg.cond[2, ,];

delta.pct.Fcond.rho  <- boot.pct.Fcond.rho - replicate(n.boot,pct.Fcond.rho(seq(l.thres,u.thres,by.thres)));
delta.pct.Fcond.pi   <- boot.pct.Fcond.pi - replicate(n.boot,pct.Fcond.pi(seq(l.thres,u.thres,by.thres)));
delta.pct.Fcond.beta <- boot.pct.Fcond.beta - replicate(n.boot,pct.Fcond.beta(seq(l.thres,u.thres,by.thres)));
delta.pct.Fcond.Fz   <- boot.pct.Fcond.Fz - replicate(n.boot,pct.Fcond.Fz(seq(l.thres,u.thres,by.thres)));
delta.pct.Fcond.selection   <- boot.pct.Fcond.selection - replicate(n.boot,pct.Fcond.selection(seq(l.thres,u.thres,by.thres)));
se.pct.Fcond.rho  <- apply(delta.pct.Fcond.rho,1,sd,na.rm=T);
se.pct.Fcond.pi   <- apply(delta.pct.Fcond.pi,1,sd,na.rm=T);
se.pct.Fcond.beta <- apply(delta.pct.Fcond.beta,1,sd,na.rm=T);
se.pct.Fcond.Fz   <- apply(delta.pct.Fcond.Fz,1,sd,na.rm=T);
se.pct.Fcond.selection   <- apply(delta.pct.Fcond.selection,1,sd,na.rm=T);
t.pct.Fcond.rho  <- apply(abs(delta.pct.Fcond.rho)/se.pct.Fcond.rho, 2, max, na.rm=T);
t.pct.Fcond.pi   <- apply(abs(delta.pct.Fcond.pi)/se.pct.Fcond.pi, 2, max, na.rm=T);
t.pct.Fcond.beta <- apply(abs(delta.pct.Fcond.beta)/se.pct.Fcond.beta, 2, max, na.rm=T);
t.pct.Fcond.Fz   <- apply(abs(delta.pct.Fcond.Fz)/se.pct.Fcond.Fz, 2, max, na.rm=T);
t.pct.Fcond.selection   <- apply(abs(delta.pct.Fcond.selection)/se.pct.Fcond.selection, 2, max, na.rm=T);

t.pct.Fcond.all <- apply(cbind(t.pct.Fcond.rho, t.pct.Fcond.pi, t.pct.Fcond.beta, t.pct.Fcond.Fz), 1, max, na.rm=T);
t.pct.Fcond.all.2 <- apply(cbind(t.pct.Fcond.selection, t.pct.Fcond.beta, t.pct.Fcond.Fz), 1, max, na.rm=T);
t.cond.joint <- apply(t.cond, c(1,3), max, na.rm=T); # n.thres*n.boot;


delta.pct.DFcond.rho  <- boot.pct.DFcond.rho - replicate(n.boot,pct.DFcond.rho);
delta.pct.DFcond.pi   <- boot.pct.DFcond.pi - replicate(n.boot,pct.DFcond.pi);
delta.pct.DFcond.beta <- boot.pct.DFcond.beta - replicate(n.boot,pct.DFcond.beta);
delta.pct.DFcond.Fz   <- boot.pct.DFcond.Fz - replicate(n.boot,pct.DFcond.Fz);
delta.pct.DFcond.selection   <- boot.pct.DFcond.selection - replicate(n.boot,pct.DFcond.selection);
se.pct.DFcond.rho  <- apply(delta.pct.DFcond.rho,1,sd,na.rm=T);
se.pct.DFcond.pi   <- apply(delta.pct.DFcond.pi,1,sd,na.rm=T);
se.pct.DFcond.beta <- apply(delta.pct.DFcond.beta,1,sd,na.rm=T);
se.pct.DFcond.Fz   <- apply(delta.pct.DFcond.Fz,1,sd,na.rm=T);
se.pct.DFcond.selection   <- apply(delta.pct.DFcond.selection,1,sd,na.rm=T);
t.pct.DFcond.rho  <- apply(abs(delta.pct.DFcond.rho)/se.pct.DFcond.rho, 2, max, na.rm=T);
t.pct.DFcond.pi   <- apply(abs(delta.pct.DFcond.pi)/se.pct.DFcond.pi, 2, max, na.rm=T);
t.pct.DFcond.beta <- apply(abs(delta.pct.DFcond.beta)/se.pct.DFcond.beta, 2, max, na.rm=T);
t.pct.DFcond.Fz   <- apply(abs(delta.pct.DFcond.Fz)/se.pct.DFcond.Fz, 2, max, na.rm=T);
t.pct.DFcond.selection   <- apply(abs(delta.pct.DFcond.selection)/se.pct.DFcond.selection, 2, max, na.rm=T);

t.pct.DFcond.all <- apply(cbind(t.pct.DFcond.rho, t.pct.DFcond.pi, t.pct.DFcond.beta, t.pct.DFcond.Fz), 1, max, na.rm=T);
t.pct.DFcond.all.2 <- apply(cbind(t.pct.DFcond.selection, t.pct.DFcond.beta, t.pct.DFcond.Fz), 1, max, na.rm=T);

delta.QFcond <- boot.QFcond - replicate(n.boot, cbind(QF2.cond(seq(l.thres,u.thres,by.thres)), 
                                                      QF2.cond.c1rho(seq(l.thres,u.thres,by.thres)), 
                                                      QF2.cond.c1rhopi(seq(l.thres,u.thres,by.thres)), 
                                                      QF2.cond.c1rhopibeta(seq(l.thres,u.thres,by.thres)), 
                                                      QF1.cond(seq(l.thres,u.thres,by.thres))));
se.QFcond <- apply(delta.QFcond, c(1,2), sd, na.rm=T);
t.QFcond.l <- lapply(1:n.boot, function(b){abs(delta.QFcond[,,b])/se.QFcond}); # n.thres*5*n.boot;
t.QFcond <- do.call(abind, list(t.QFcond.l,along=3));
t.QFcond.joint <- apply(t.QFcond, c(1,3), max, na.rm=T); # n.thres*n.boot;



crt.cond <- apply(apply(t.cond[qn.lower:qn.upper,,],c(2:3),max,na.rm=T),1,quantile,1-alpha,na.rm=T); # 5*1
crt.avg.cond <- apply(t.avg.cond, 1, quantile, 1-alpha, na.rm=T);
crt.pct.Fcond.rho  <- quantile(t.pct.Fcond.rho, 1-alpha, na.rm=T);
crt.pct.Fcond.pi   <- quantile(t.pct.Fcond.pi, 1-alpha, na.rm=T);
crt.pct.Fcond.beta <- quantile(t.pct.Fcond.beta, 1-alpha, na.rm=T);
crt.pct.Fcond.Fz   <- quantile(t.pct.Fcond.Fz, 1-alpha, na.rm=T);
crt.pct.Fcond.selection <- quantile(t.pct.Fcond.selection, 1-alpha, na.rm=T);

crt.pct.Fcond.all <- quantile(t.pct.Fcond.all, 1-alpha, na.rm=T);
crt.pct.Fcond.all.2 <- quantile(t.pct.Fcond.all.2, 1-alpha, na.rm=T);
crt.cond.joint <- quantile(apply(t.cond.joint[qn.lower:qn.upper,],2,max,na.rm=T), 1-alpha, na.rm=T);
crt.avg.cond.joint <- quantile(apply(t.avg.cond, 2, max, na.rm=T), 1-alpha, na.rm=T);

crt.pct.DFcond.rho  <- quantile(t.pct.DFcond.rho, 1-alpha, na.rm=T);
crt.pct.DFcond.pi   <- quantile(t.pct.DFcond.pi, 1-alpha, na.rm=T);
crt.pct.DFcond.beta <- quantile(t.pct.DFcond.beta, 1-alpha, na.rm=T);
crt.pct.DFcond.Fz   <- quantile(t.pct.DFcond.Fz, 1-alpha, na.rm=T);
crt.pct.DFcond.selection <- quantile(t.pct.DFcond.selection, 1-alpha, na.rm=T);

crt.pct.DFcond.all <- quantile(t.pct.DFcond.all, 1-alpha, na.rm=T);
crt.pct.DFcond.all.2 <- quantile(t.pct.DFcond.all.2, 1-alpha, na.rm=T);

crt.QFcond.joint <- quantile(apply(t.QFcond.joint[qn.lower:qn.upper,],2,max,na.rm=T), 1-alpha, na.rm=T);

Fcond <- cbind(F2.cond, F2.cond.c1rho, F2.cond.c1rhopi, F2.cond.c1rhopibeta, F1.cond);
e.avg.Fcond <- c(e.avg.F2.cond, e.avg.F2.cond.c1rho, e.avg.F2.cond.c1rhopi, e.avg.F2.cond.c1rhopibeta, e.avg.F1.cond);


ubound.cond <- apply(Fcond + t(crt.cond*sqrt(var.cond)), 2, sort);
ubound.cond[ubound.cond > 1] <- 1;
ubound.cond[ubound.cond < 0] <- 0;
lbound.cond <- apply(Fcond - t(crt.cond*sqrt(var.cond)), 2, sort);
lbound.cond[lbound.cond > 1] <- 1;
lbound.cond[lbound.cond < 0] <- 0;

ubound.cond.joint <- apply(Fcond + t(crt.cond.joint*sqrt(var.cond)), 2, sort);
ubound.cond.joint[ubound.cond.joint > 1] <- 1;
ubound.cond.joint[ubound.cond.joint < 0] <- 0;
lbound.cond.joint <- apply(Fcond - t(crt.cond.joint*sqrt(var.cond)), 2, sort);
lbound.cond.joint[lbound.cond.joint > 1] <- 1;
lbound.cond.joint[lbound.cond.joint < 0] <- 0;


ubound.avg.cond <- e.avg.Fcond + crt.avg.cond*sqrt(var.avg.cond);
lbound.avg.cond <- e.avg.Fcond - crt.avg.cond*sqrt(var.avg.cond);
ubound.avg.cond.joint <- e.avg.Fcond + crt.avg.cond.joint*sqrt(var.avg.cond);
lbound.avg.cond.joint <- e.avg.Fcond - crt.avg.cond.joint*sqrt(var.avg.cond);

ubound.pct.cond.rho  <- function(x){pct.Fcond.rho(x) + crt.pct.Fcond.all*
                                      approxfun(se.pct.Fcond.rho, x=seq(l.thres,u.thres,by.thres), rule=2)(x)};
lbound.pct.cond.rho  <- function(x){pct.Fcond.rho(x) - crt.pct.Fcond.all*
                                      approxfun(se.pct.Fcond.rho, x=seq(l.thres,u.thres,by.thres), rule=2)(x)};
ubound.pct.cond.pi   <- function(x){pct.Fcond.pi(x) + crt.pct.Fcond.all*
                                      approxfun(se.pct.Fcond.pi, x=seq(l.thres,u.thres,by.thres), rule=2)(x)};
lbound.pct.cond.pi   <- function(x){pct.Fcond.pi(x) - crt.pct.Fcond.all*
                                      approxfun(se.pct.Fcond.pi, x=seq(l.thres,u.thres,by.thres), rule=2)(x)};
ubound.pct.cond.beta <- function(x){pct.Fcond.beta(x) + crt.pct.Fcond.all*
                                      approxfun(se.pct.Fcond.beta, x=seq(l.thres,u.thres,by.thres), rule=2)(x)};
lbound.pct.cond.beta <- function(x){pct.Fcond.beta(x) - crt.pct.Fcond.all*
                                      approxfun(se.pct.Fcond.beta, x=seq(l.thres,u.thres,by.thres), rule=2)(x)};
ubound.pct.cond.Fz   <- function(x){pct.Fcond.Fz(x) + crt.pct.Fcond.all*
                                      approxfun(se.pct.Fcond.Fz, x=seq(l.thres,u.thres,by.thres), rule=2)(x)};
lbound.pct.cond.Fz   <- function(x){pct.Fcond.Fz(x) - crt.pct.Fcond.all*
                                      approxfun(se.pct.Fcond.Fz, x=seq(l.thres,u.thres,by.thres), rule=2)(x)};
ubound.pct.cond.selection.2  <- function(x){pct.Fcond.selection(x) + crt.pct.Fcond.all.2*
                                            approxfun(se.pct.Fcond.selection, x=seq(l.thres,u.thres,by.thres), rule=2)(x)};
lbound.pct.cond.selection.2  <- function(x){pct.Fcond.selection(x) - crt.pct.Fcond.all.2*
                                            approxfun(se.pct.Fcond.selection, x=seq(l.thres,u.thres,by.thres), rule=2)(x)};
ubound.pct.cond.beta.2 <- function(x){pct.Fcond.beta(x) + crt.pct.Fcond.all.2*
                                      approxfun(se.pct.Fcond.beta, x=seq(l.thres,u.thres,by.thres), rule=2)(x)};
lbound.pct.cond.beta.2 <- function(x){pct.Fcond.beta(x) - crt.pct.Fcond.all.2*
                                      approxfun(se.pct.Fcond.beta, x=seq(l.thres,u.thres,by.thres), rule=2)(x)};
ubound.pct.cond.Fz.2   <- function(x){pct.Fcond.Fz(x) + crt.pct.Fcond.all.2*
                                      approxfun(se.pct.Fcond.Fz, x=seq(l.thres,u.thres,by.thres), rule=2)(x)};
lbound.pct.cond.Fz.2   <- function(x){pct.Fcond.Fz(x) - crt.pct.Fcond.all.2*
                                      approxfun(se.pct.Fcond.Fz, x=seq(l.thres,u.thres,by.thres), rule=2)(x)};

ubound.pct.DFcond.rho  <- pct.DFcond.rho + crt.pct.DFcond.all*se.pct.DFcond.rho;
lbound.pct.DFcond.rho  <- pct.DFcond.rho - crt.pct.DFcond.all*se.pct.DFcond.rho;
ubound.pct.DFcond.pi   <- pct.DFcond.pi + crt.pct.DFcond.all*se.pct.DFcond.pi;
lbound.pct.DFcond.pi   <- pct.DFcond.pi - crt.pct.DFcond.all*se.pct.DFcond.pi;
ubound.pct.DFcond.beta <- pct.DFcond.beta + crt.pct.DFcond.all*se.pct.DFcond.beta;
lbound.pct.DFcond.beta <- pct.DFcond.beta - crt.pct.DFcond.all*se.pct.DFcond.beta;
ubound.pct.DFcond.Fz   <- pct.DFcond.Fz + crt.pct.DFcond.all*se.pct.DFcond.Fz;
lbound.pct.DFcond.Fz   <- pct.DFcond.Fz - crt.pct.DFcond.all*se.pct.DFcond.Fz;

ubound.pct.DFcond.selection.2 <- pct.DFcond.selection + crt.pct.DFcond.all.2*se.pct.DFcond.selection;
lbound.pct.DFcond.selection.2 <- pct.DFcond.selection - crt.pct.DFcond.all.2*se.pct.DFcond.selection;
ubound.pct.DFcond.beta.2 <- pct.DFcond.beta + crt.pct.DFcond.all.2*se.pct.DFcond.beta;
lbound.pct.DFcond.beta.2 <- pct.DFcond.beta - crt.pct.DFcond.all.2*se.pct.DFcond.beta;
ubound.pct.DFcond.Fz.2   <- pct.DFcond.Fz + crt.pct.DFcond.all.2*se.pct.DFcond.Fz;
lbound.pct.DFcond.Fz.2   <- pct.DFcond.Fz - crt.pct.DFcond.all.2*se.pct.DFcond.Fz;

ubound.QFcond.joint <- cbind(QF2.cond(seq(l.thres,u.thres,by.thres)),
                             QF2.cond.c1rho(seq(l.thres,u.thres,by.thres)),
                             QF2.cond.c1rhopi(seq(l.thres,u.thres,by.thres)),
                             QF2.cond.c1rhopibeta(seq(l.thres,u.thres,by.thres)),
                             QF1.cond(seq(l.thres,u.thres,by.thres))) + crt.QFcond.joint*se.QFcond;
lbound.QFcond.joint <- cbind(QF2.cond(seq(l.thres,u.thres,by.thres)),
                             QF2.cond.c1rho(seq(l.thres,u.thres,by.thres)),
                             QF2.cond.c1rhopi(seq(l.thres,u.thres,by.thres)),
                             QF2.cond.c1rhopibeta(seq(l.thres,u.thres,by.thres)),
                             QF1.cond(seq(l.thres,u.thres,by.thres))) - crt.QFcond.joint*se.QFcond;

ubound.QF.cond.joint <- apply(ubound.cond.joint, 2, approxfun, y=ys, rule=2);
lbound.QF.cond.joint <- apply(lbound.cond.joint, 2, approxfun, y=ys, rule=2);


################
# Save Results #
################
save.image(paste("Results",str.sample,str.rho,str.trim,str.cb,str.robust,str.boot,".RData",sep=""));


